#include "ActiveRegion.h"
#include "Matrix.h"
#include "Global.h"
#include <string>
#include <fstream>
#include <sstream>
#include <iostream>
#include <cmath>

using namespace std;

ActiveRegion::ActiveRegion(char * filename)
{
  int i;
  ifstream inputFile;
  string line;
 
  inputFile.open(filename, ios::in);

  grid_size = -1;
  while ( inputFile.good() )
  {
    getline(inputFile, line);
    grid_size++;
  }
  
  x_grid = new double[grid_size];
  V = new double[grid_size];

  inputFile.clear();
  inputFile.seekg(0, ios::beg);

  //assign grid and V
  for (i = 0; i < grid_size; i++)
  {
    getline(inputFile, line, ',');
    istringstream( line ) >> x_grid[i];
    getline(inputFile, line);
    istringstream(line) >> V[i];
  }

  inputFile.close();

  delta_x = x_grid[1] - x_grid[0];

}


ActiveRegion::~ActiveRegion()
{
  delete [] x_grid;
  delete [] V;

}

//double operator/(const double left[], const double right[]);


void ActiveRegion::MakeHamiltonian()
{
	double A = (hbar*hbar)/(2*mass*echarge*(delta_x*delta_x));
	Matrix Hamilt(grid_size);

	//MAKE HAMILTONIAN
	for(int i = 0; i < grid_size; i++)
	{
		Hamilt[i][i] = (2*A + V[i]);
			if(i == 0)
				Hamilt[1][0] = -1*A;
			else if(i == (grid_size - 1))
				Hamilt[i-1][i] = -1*A;
			else
			{
				Hamilt[i-1][i] = -1*A;
				Hamilt[i+1][i] = -1*A;
			}
	}
	
	//Find 0th Eigenfunctions
	Fork0 = new double [grid_size];
	for(int i = 0; i < grid_size; i++)
		Fork0[i] = 1;


	for(int x = 0; x < 5000; x++)
	{
		Fork0 = Hamilt*Fork0;
	
		double sum = 0;
		for(int n = 0; n <grid_size; n++)	//NORMALIZE
		{
			sum = sum + Fork0[n]*Fork0[n]*delta_x;
		}

		for(int n = 0; n <grid_size; n++)
		{
			Fork0[n] = Fork0[n] /sqrt(sum);
		}
		cout<<x;
		if(x == 4999)
			cout << endl << "SUM: " << sum;
	}
	cout << endl;
	for(int n = 0; n <grid_size; n++)
		{
			//cout << Fork0[n] << endl;
		}
	//search for greatest
	double greatest = 0;
	int index;

	for(int i = 0; i < grid_size; i++)
	{
		if(Fork0[i] > greatest)
		{
			greatest = Fork0[i];
			index = i;
		}
	}

	cout << endl << "INDEX: " << index << "	FORK0[index]: " << greatest << endl;

	
	/*double* Fork1 = new double [grid_size];
	for(int i = 0; i < grid_size; i++)
		Fork1[i] = delta_x;

	for(int x = 0; x < 5000; x++)
	{
		Fork0 = Hamilt*Fork0;
	
		double sum = 0;
		for(int n = 0; n <grid_size; n++)	//NORMALIZE
		{
			sum = sum + Fork0[n]*Fork0[n]*delta_x;
		}

		for(int n = 0; n <grid_size; n++)
		{
			Fork0[n] = Fork0[n] /sqrt(sum);
		}
		cout<<x;
		if(x == 4999)
			cout << endl << "SUM: " << sum;
	}*/
}